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THE  POSTPROCESSING  TECHNIQUE  IN  THE  FINITE  ELEMENT  METHOD. 

THE  THEORY  AND  EXPERIENCE 

I .  Babuska^ 

University  of  Maryland,  College  Park,  Maryland,  U.S.A. 

2 

K.  Izadpanah 
B.  Szabo^ 

Washington  University,  St.  Louis,  Missouri,  U.S.A. 


The  paper  addresses  the  h,  p,  and  h-p  versions 
of  the  finite  element  method  in  connection  with  a 
postprocessing  technique  for  extracting  the  values 
of  a  functional.  This  technique  combines  the 
finite  element  method  with  the  analytical  ideas  of 
the  theory  of  partial  differential  equations  of 
elliptic  type. 


1 .  INTRODUCTION 

Finite  element  computations  in  structural  mechanics  usually 
have  two  purposes:  (1)  to  determine  the  stress  and 
displacement  fields  and  (2)  to  determine  the  values  of  certain 
functionals  defined  on  displacement  fields  as,  for  example,  the 
stress  intensity  factors,  stresses  at  specific  points, 
reactions,  etc.  Computations  of  these  values  involve  the 
finite  element  solution.  For  example,  the  stress  components 
are  often  computed  at  the  Gauss  points  of  the  elements  and  the 
stresses  at  any  other  points  are  then  computed  by  the 
interpolation  technique,  the  stress  intensity  factors  is 
determined  by  the  J-integral  or  curve  fitting  technique, 
etc.  We  shall  refer  to  these  operations  as  postprocessing. 

Usually  the  values  of  these  functionals  are  needed  to  be  known 
with  higher  accuracy  and  reliability  than  the  displacement  or 
stress  field  itself. 


^Partially  supported  by  the  Office  of  Naval  Research  under 
grant  number  N0001 4- 77-C-0623 . 

^Partially  supported  by  the  Office  of  Naval  Research  under 
grant  number  N0001 4-81 -K-0625 . 
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Assuming  that  we  have  the  finite  element  solution  and  wish  to 
determine  certain  functional  values  the  following  questions 
arise: 

1 )  What  should  the  relationship  be  between  the  computa¬ 
tional  effort  spent  on  the  finite  element  solution  and  the 
effort  spent  on  postprocessing:  Is  it  better  to  use  a  very- 
simple  and  inexpensive  postprocessing  technique  as  for  example 
direct  evaluation  of  the  stresses  from  the  finite  element 
solution  in  the  desired  points  or  should  one  select  a  more 
expensive  technique.  Of  course  we  have  to  relate  the  answer  to 
the  achieved  accuracy  and  to  the  reliability  and  robustness  of 
the  postprocessing  procedures  under  consideration. 

2)  Given  a  finite  element  solution,  what  is  the  largest 
accuracy  of  the  functional  values  one  can  achieve  by  the 
postprocessing  technique.  In  other  words,  what  is  the  maximal 
information  contained  in  the  finite  element  solution  which 
could  be  used  for  the  extraction  of  the  desired  value. 

3)  How  do  the  various  versions  of  the  finite  element 
method,  i.e.,  the  h-version,  the  p-version  and  the  h-p 
version  bear  on  the  importance  of  proper  selection  of  the 
postprocessing  techniques. 

These  questions  are  discussed  in  some  details. 


2.  THE  EXTENSION  OPERATORS.  THE  h,  p  AND  h-p  VERSIONS  OF 
THE  FINITE  ELEMENT  METHOD 

There  are  three  versions  of  the  finite  element  methods  based  on 
the  common  variational  (energy)  principle.  They  are  charac¬ 
terized  by  the  systematic  selection  (extension)  of  the  finite 
element  spaces  leading  to  the  convergence  of  the  finite  element 
solutions  to  the  exact  one. 

The  h-version  is  the  classical  and  most  commonly  used  method  of 
extension:  the  polynomial  degree  of  elements  p  is  fixed  and 

mesh  refinement  is  used  for  controlling  the  error  of  approxi¬ 
mation  (h  refers  to  the  size  of  the  element).  Typically  the 
polynomial  degree  of  elements  is  low,  usually  p  =  1  or  2. 
Proper  selection  of  the  mesh  and  its  refinement  strongly 
influences  the  error  and  its  dependence  on  the  computational 
effort. 

In  the  p-version  the  mesh  is  fixed  and  the  polynomial  degree  of 
elements  is  increased  either  uniformly  or  selectively  over  the 
mesh. 

The  h-p  version  combines  the  h  and  p-versions,  i.e.,  error 
reduction  is  achieved  by  a  proper  mesh  refinement  and  con¬ 
current  changes  in  the  distribution  of  the  polynomial  degree  of 
elements. 


The  performance  of  the  various  extension  operators  can  be 
compared  from  various  points  of  view,  the  most  important  of 
which  are  human  and  computer  resource  requirements  in  relation 
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to  the  desired  level  of  precision.  Such  relationships  are 
difficult  to  quantify  and  are  subject  due  to  various  factors, 
therefore  the  performance  of  the  extension  operators  is  usually 
related  to  the  number  of  degrees  of  freedom  N.  Of  course 
evaluation  of  an  extension  operator  would  not  be  meaningful 
without  considering  the  goals  of  computation.  For  example,  if 
only  stress  intensity  factors  are  desired,  then  the  accuracy  of 
the  computed  displacements,  reactions  or  stresses  are  not  of 
importance.  In  many  cases  the  computation  has  multiple  goals. 


3.  THE  MODEL  PROBLEM 


In  order  to  illustrate  the 


Figure  1 

The  model  problem 


essential  properties  of  finite 
element  solution  and 
extraction  techniques,  we 
have  selected  a  model 
problem  which  represents 
some  of  key  features  of  a 
large  class  of  engineering 
problems.  Specifically  let 
us  consider  the  plane 
strain  problem  of  two- 
dimensional  elasticity 
(homogeneous  isotropic 
material)  with  E  and  v 
representing  the  modulus  of 
elasticity  and  Poisson 
ratio  respectively  (E  >  0, 
0  <  v  <  .5).  The  domain 
D,  a  square  panel  with  a 
crack,  is  shown  in  Fig.  1  . 


We  shall  be  concerned  here  with  problems  in  which  only 
tractions  are  prescribed  at  the  boundary  (i.e.,  first  boundary 
value  problem  of  elasticity). 

We  denote  the  displacement  vector  function  by  u  =  {u. ,u?} 
and  the  corresponding  stress  tensor  by  “  " 
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The  solution  u  satisfies  the  Navier-Larae  equations.  It  is 
possible  to  express  the  solution  through  two  holomorphic 
functions  $(z),  ♦(z)  using  the  theory  of  Muskhelishvili  [ij. 


2u(u^  +  iu2)  =  «$(z)  -  zTTzT  -  TTzJ  (3-2) 

where 

z  =  Xi  +  ix2,  y  =  ,  <  =  3  -  4v  (3.3) 


and  z  =  x1  -  iX2»  resp.  $’(z)  mean  conjugate  values  to  z 
and  ♦ ' ( z ) . 

The  components  of  the  stress  tensor  are  expressed  by  Kolosov- 
Muskhelishvili  formulae 


T11  +  t22  =  +  TTzT)  =  4  Re  4>  *  ( z ) 

=  2(*(z)  +  •'(  z) )  (3.4) 


t22  “  T11  +  2iT12  =  2tz*"^z)  +  =  2[ z* ' ( z)  +  T(z)J 


where 


*(z)  =  ♦'(z),  *(z)  =  t|> '  ( z ) 


(3.5) 


(3.6) 


and  Re  ♦’(z)  is  the  real  part  of  ♦'(z). 

The  correspondence'  between  the  displacements  (and  the  stress) 
field  and  the  functions  $  and  is  one  to  one  up  to  the 

constants  y  and  y'  in  $  and  \|»,  respectively,  satisfying 
the  relation  y  -  y'  =0. 


In  our  model 


problem  we  consider  the 
*(z)  =  (1+i)z_ 

n(z)  =  «(z) 

n(z)  = 


following  (exact) 


solution 

(3.7) 

(3.8) 

(3.9) 


M\ 


*(z)  +  z*'(z)  +  *p(z) 
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where  *(z)  =  *  ( z )' ,  *’(z)  =  *'(z7,  y(z)  =  H'(z). 

fl( z)  is  a  holomorphic  function  on  D.  ^unction  z~  ^2  is  to 

be  understood  as  the  principal  branch  of  z”  ^2  on  T). 
Function  ?(z)  is  uniquely  defined  by  (3*9)  and  (3*7)  (3.8). 
The  tractions  on  the  boundary  of  D  are  defined  by  (3*4) 

(3*5) •  It  can  be  readily  verified  that  the  two  edges  of  the 
crack  are  traction  free. 


We  will  now  discuss  the  finite  element  solution  and  the 
postprocessing  technique  if  the  tractions  are  prescribed  on  the 
boundary  of  D  so  that  the  exact  solution  to  the  problem  is 
given  by  (3.7)  —  (3-9)  •  Specif icaly  we  now  consider  the  case  E 
=  1 ,  v  =  3*  The  strain  energy  of  the  exact  solution  is:  W  = 
42.16491240. 


4.  THE  FINITE  ELEMENT  SOLUTION 

We  have  solved  the  model  problem  by  the  h  and  p-versions  of 
the  finite  element  method.  The  p-version  of  the  finite 


Figure  2 

The  meshes  for  the  p-version,  A:  Mesh  1 ,  B:  Mesh  2 


element  method  was  implemented  in  the  experimental  computer 
program  COMET-X  developed  at  the  Center  for  Computational 
Mechanics  of  Washington  University  in  St.  Louis  [2j.  The  two 
meshes  shown  in  Fig.  2A,B  were  used.  The  polynomial  degrees 
were  the  same  for  all  elements  and  ranged  from  1  to  8.  The 
shape  functions  on  trapezoidal  elements  of  mesh  2  were 
constructed  by  blending  function  technique. 


The  h-version  solution  was  obtained  by  means  of  the  computer 
program  FEARS  developed  at  the  University  of  Maryland  [3,  4j. 


cm 


m 


vl 


PEARS  uses  quadrilateral  elements  of  degree  one.  The  program 
is  adaptive  and  produces  a  sequence  of  nearly  optimal  meshes. 
See  [3J  [4]  [5]  [6]  [7].  The  mesh  from  this  sequence  with  319 
elements  and  number  of  degrees  of  freedom  N  =  617  is  shown  in 
Pig.  3. 


Figure  3 

The  mesh  constructed  by  the  adaptive  program  PEARS 


5.  ERROR  OP  THE  FINITE  ELEMENT  SOLUTION  MEASURED  IN  ENERGY 
NORM 


We  denote  the  exact  solution  by  uQ  and  the  finite  element 
solution  by  upj>*  The  error  of  the  finite  element  solution  is 

denoted  by  s> 


§  *  Ho  "  Hpe* 


We  measure  the  magnitude  of  the  error  by  the  energy  norm  i*igt 


lei,  -  [W(e)]V2  . 


(5.1  ) 
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This  measure  is  equivalent  to  measuring  the  error  in  the  stress 
components  by  integrals  of  their  squares  ( the  L2  norm).  In  our 
case  when  tractions  are  specified  at  the  boundary 

W(upE)  <  W(uQ)  (5.2) 

and 

ieiE  =  [W(uQ)  -  W(upE)]^2  .  (5. "5) 

The  extension  operators  under  consideration  monotonically 
increase  the  finite  element  spaces  either  by  increasing  the 
degree  of  elements  or  refining  the  mesh.  Therefore  the  energy 
norm  of  the  error  monotonically  decreases.  We  can  write 

ieiE  i  C(N)N"W  (5-4) 

and  expect  that  for  properly  chosen  p  the  function  C(N)  is 
nearly  constant  especially  for  larger  N.  The  number  p  >  0 
is  the  rate  of  convergence  of  the  error  measured  in  the  energy 
norm. 

It  is  possible  to  estimate  the  value  of  p.  In  our  case  the 
rate  p  is  governed  by  the  strength  of  the  singularity  of  the 
solution.  It  can  be  shown  that  for  the  p-version  [8];[9] 

I  el E  <  C(e)N"(  1/2  "e)  (5.5) 

with  e  >  0  arbitrarily  small  and  C  independent  of  N. 

The  h-version  using  the  uniform  mesh  yields  the  estimate 

1  el E  <  CN~  %  (5.6) 

with  the  rate  independent  of  the  degree  of  elements.  The 
optimal  refinement  of  the  mesh  leads  to  the  estimate 

» el  E  <  CN“p/2  (5.7) 


(PEARS  uses  p  =  1)  where  the  rate  is  independent  of  the 
strength  of  the  singularity. 


The  h-p  version  with  optimal  mesh  and  p-distribution  leads  to 
the  estimate 


'2'e 


<  Ce 


where  8  =  1  /3  independently  of  the  strength  of  the 
singularity  and  y  >  0. 

The  relative  error  in  the  energy  norm  defined  as 


8 


'-'e.r 


■Ho'e 


(5.8) 


has  been  plotted  in  Pig.  4  on  log— log  scale  for  the  p— version 
(mesh  1,2),  for  the  h-version  with 
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Figure  4 

Relative  error  in  the  energy  norm  vs  degrees  of  freedom 

(1)  h-version,  uniform  mesh,  (2)  h-version,  adaptively 
constructed  mesh,  (3)  p-version  Mesh  1,  (4)  p-version 
Mesh  2 


adaptively  constructed  mesh  and  for  the  h-version  with  uniform 
mesh.  The  polynomial  degree  of  elements  is  also  shown  in  the 
figure.  The  shown  slopes  are  the  theoretical  slopes  of  the 
rate  of  convergence  [m  =  V2  and  V4  ].  It  is  seen  that  the 
observed  rate  of  convergence  closely  agrees  with  (5. 5) (5. 7)  • 
Prom  (5.4)  we  can  compute  C(N)  for  the  p-version.  The 
results  are  given  in  Table  1 . 

Tables  2  and  3  show  analogous  results  for  the  h-version.  The 
comparison  between  Tables  1-3  shows  that  for  5%  accuracy  we 
need  N  *  1770  when  using  p-version  Mesh  2,  N  =  2290  for 
h-version  with  adaptively  refined  mesh  and  N  *  146000  for  h- 
version  with  uniform  mesh. 


TABLE  1 


ca.ii 


Relationship  between  (e<_  D  and  N  for  the 

ft ,  K 

p-version,  Mesh  2  [p  =  '/?  J 


C  ( M ) / lu-  I  _ 


p 

N 

1 

35 

2 

95 

3 

135 

4 

239 

5 

347 

6 

479 

7 

635 

8 

815 

Mill 


1  .816 
1  .997 
2.059 
2.061 
2.079 
2.088 
2.099 


TABLE  2 


Relationship  between  Help,  R  and  N  for  the  h-version 
with  adaptively  constructed  mesh  [p  =  Vp  J 


C(N)/iu0«e 


ittVI 


2.665 

2.562 

2.501 

2.366 

2.394 


TABLE  3 


Relationship  between  lei™  n  and  N  for  the  h 
version  with  uniform  mesh  [p  =  '4  J 


C(R)/l  Uq  i  pi 


36 .02# 
21.01% 
19-81* 


.967 

.974 

.977 


6.  COMPUTATION  OP  THE  STRESSES 

The  finite  element  method  provides  the  solution  u-ppi  which 
converges  to  the  exact  solution  in  the  energy  norm.  We  have 
seen  that  the  error  measured  in  this  norm  decreases  raonoton- 
ically  and  in  a  very  orderly  way.  We  now  examine  the  pointwise 
error  an  stresses  for  the  h  and  p-versions.  We  denote  the 
error  in  the  stress  components  as 


?  V 


10 


[Oj,  v  [ PE  J /  v 

•ij  (X,,X2)  -  Tj.  (X,,X2) 


(6.1  ) 


and  the  relative  error  by 


(x,.V  - 

|tf°J(x, ,x2) 


r  o  l  r  I 

where  and  tV^.  J  are  respectively  the  stress  compo¬ 

nents  corresponding  to  the  exact  and  finite  element  solutions. 
We  will  compute  the  stresses  directly  from  the  derivatives  of 
u^g  and  the  stress-strain  law.  Pig.  5  shows  the  relative  erior 

e.^  in  at  the  point  (.0,.1)  computed  by  the  p-version. 
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Figure  5 

The  relative  error  of  e^j  computed  by  the  p-version 
a)  Mesh  1,  b)  Mesh  2.  (1)  e^  ,  (2)  e^,  (3)  e*2 


Pig.  6  shows  isometric  drawings  of  the  error  in  f°r 

various  p  values  for  meshes  1  and  2.  The  error  values  were 
computed  on  a  uniform  grid  with  the  grid  points  (ih,jh)  h  = 

.1 ,  i,j  =  -  10,  10.  At  points  other  than  the  grid  points,  the 
values  were  computed  by  linear  interpolation. 

In  the  case  of  the  h-version,  the  error  is  discontinuous  at 
the  boundary  of  every  element.  Therefore  we  compute  the 
stresses  in  the  center  of  every  element  where  increased 
accuracy  can  be  expected. 

In  Pig.  7  we  show  the  level-lines  of  the  error  in  t ^  (using 

the  mesh  shown  in  Pig.  3)  in  the  upper  right  quarter  of  the 
domain  D.  The  local  maxima  and  minima  are  shown  also  in  the 
figure.  The  error  is  large  in  the  neighborhood  of  the  tip  of 
the  crack.  The  level-lines  and  the  local  maxima  and  minima 
depend  on  the  interpolation  technique  used.  We  see  in  contrast 
to  the  p-version  that  the  oscillatory  behaviour  of  the  error 
is  not  so  strong  here;  nevertheless,  it  has  to  be  emphasized 
that  when  the  stresses  are  computed  everywhere  directly  from 
the  displacements,  strong  oscillatory  behaviour  will  appear  in 
every  element. 

The  center  of  the  elements  are  changing  with  the  mesh.  To 
show  the  convergence  of  the  stresses,  we  selected  for  the  Table 
4  the  center  points  which  are  closest  to  the  tip  of  the  crack 
(in  the  first  quarter  of  D).  The  table  shows  the  error  in  < 
and  the  magnitude  of  the  exact  values  of  the  stress. 


Pigure  7 

The  level  lines  of  the  error  e^  of 
quarter  of  T)  computed  by  the  h-version 


in  the  upper 


TABLE  4 


The  relative  error  of  the  stresses  in  the  neighborhood  of  the 
origin. 


No.  of 
elements 

N 

Coordinates 

leR  I 
e11  1 

1  e22 

leE  1 
e12 

*1 

*2 

,T11  1 

t22  1 

|T[0]| 
t12  1 

■ 

16 

51 

.25 

.25 

35.54# 

5.037 

19.90# 

3.751 

20.70# 

1.553 

Adaptive  Mesh 

43 

143 

.125 

.125 

33.95# 

7.124 

10.84# 

5.304 

15.14# 

2.197 

106 

221 

•3125(-1) 

«3125(-1 ) 

31 .01# 
14.251 

6.09# 

10.609 

12.35# 

4.394 

319 

617 

.71 25  (-2) 

.71 25  (-2) 

29-77# 

28.501 

2.40# 

21.218 

17.65# 

8.789 

16 

51 

.25 

.25 

35.54# 

5.037 

19.90# 

3.751 

20.70# 

1.553 

Uniform  Mesh 

64 

167 

.125 

.125 

36.34# 

7.124 

20.07# 

5.304 

14.52# 

2.197 

256 

591 

.625 (-2) 

.625 (-1 ) 

34.46# 

10.076 

17.16# 

7.502 

14.09# 

3-107 

To  depict  the  behaviour  in  a  fixed  point  (.25,  .25)  we  select 
the  center  points  closest  to  it.  Table  5  shows  the  results. 

If  we  desire  to  compute  the  stress  components  in  the  nodal 
point  (.25,  .25)  we  have  4  values  for  disposition  and  also 
their  average.  Table  6  we  shows  the  relative  errors.  The 
value  in  the  lines  1,  2,  3»  4  are  computed  from  the  elements 
ordered  counterclockwise  starting  with  the  upper-right  one. 

The  line  A  shows  relative  error  of  the  average  of  the  stress 
values  computed  in  the  four  elements. 

In  contrast  to  the  monotonic  and  orderly  behaviour  of  the  error 
measured  in  the  energy  norm,  the  accuracy  in  the  stresses  is 
poor  and  nonmonotonic,  although  the  stresses  are  converging  in 
integral  sense  (in  the  energy  norm)  raonotonically .  In 
addition,  the  quality  of  the  computed  stress  components  is  very 
different. 


TABLE  5 


i 

i 


The  relative  error  of  the  stresses  in  the  neighborhood  of 
(.25, .25). 


No.  of 
elements 

Coordinates 

le^  1 
e22 

mm 

■ 

N 

X1 

*2 

,Tto] 

,Tii 

leR  1 
e12 

■ 

16 

51 

.25 

.25 

35.54# 

5.037 

19-90# 

3.751 

20.70# 

1.553 

Adaptive  Mesh 

43 

143 

.375 

.375 

1.09# 

4.113 

2.10# 

3.062 

43.51# 

1 .268 

106 

221 

.1875 

.1875 

3.89# 

5.817 

8.26# 

4.331 

40.28# 

1.794 

319 

617 

.21875 

.21875 

.464# 

5.386 

.812# 

4.010 

12.12# 

1.661 

16 

51 

.25 

.25 

35.54# 

5.037 

19.90# 

3.751 

20.70# 

1.553 

Uniform  Mesh 

64 

167 

.373 

.375 

8.01# 

4.113 

7.79# 

3.062 

16.06# 

1.268 

256 

591 

.1875 

.1875 

10.20# 

5.817 

10.78# 

4.331 

6.45# 

1.794 

7.  POSTPROCESSING 

We  have  seen  that  stresses  computed  directly  from  finite 
element  solutions  are  not  accurate.  Nevertheless,  often  the 
values  of  the  stresses  is  the  main  aim  of  the  computation. 

We  will  show  now  that  by  utilizing  the  analytical  structure  of 
the  Navier-Lame  equations  it  is  possible  to  compute  stresses 
with  the  accuracy  comparable  to  the  accuracy  of  the  energy  of 
the  finite  element  solution  (which  is  the  square  of  the  error 
measured  in  the  energy  norm).  We  will  outline  the  main  idea. 
For  more,  see  [8],  [9],  [10] . 

Let  Xq  =  ( Xq , i » xo , 2 ^  (  D  and  denote  by  S(xQ,p)  the  disc  of 
radius  p  centered  in  xQ.  Further,  let  D(xQ,p)  =  D 
-S(xQ,p).  See  Fig.  8.  The  boundary  of  D(xQ,p)  is  denoted 
by  9D(Xq,p)  =  an  U  r  where  r  is  the  boundary  of  the  disk 
S(x0,p)  .  We  now  define  the  extraction  (displacement)  function 
w(xQ,x)  s  (w1#w2)  which  corresponds  to  the  functions  +  ,  J  in 


No.  of 
elements 

N 

43 

143 

106 

221 

319 

617 

64 

167 

256 

591 

eR* 

eH% 

e22  * 

R  j 
e12* 

.034 

1.60 

33.95 

.042 

4.97 

4.35 

2.41 

6.60 

3.38 

.026 

1.76 

72.22 

3.09 

3.14 

64.41 

10.99 

7.57 

11.16 

4.79 

2.31 

35.86 

12.21 

2.01 

.093 

17.27 

17.42 

106.43 

9.71 

13.01 

71 .49 

4.09 

5.47 

13.42 

1.46 

2.68 

22.96 

4.41 

.099 

55.69 

4.41 

13.43 

3.88 

6.07 

11.74 

28.85 

12.12 

10.65 

17.36 

19.37 

15.35 

13.36 

6.75 

8.68 

5.47 

5.87 

5.94 

68.41 

17.51 

12.63 

60.22 

8.11 

10.13 

12.66 

8.10 

6.62 

13.38 

5.38 

5.05 

8.15 

8.12 

13.64 

38.71 

10.84 

15.84 

33.48 

i*v  v" «  '.cso 


I  *,»  VfAtXi 


Uniform  Mesh  I  Adaptive  Mesh 


displacement  function  w  defined  by  (7.1)  through  (7.3)  by 
(3.2)  is  a  single  valued  function  and  it  is  an  admissible 
displacement  function. 


Denote  by  the  stress  tensors  associated  with  the 


Pig.  8.  The  domain  D(xn,p). 


displacement  functions  u  and  w.  Denote  the  outward  normal  to 
9ft (xq , p)  by  n.  Then  Betti’s  law  can  be  written  in  the  form 


/  (u,T^«n)ds  =  / 

3D(xq,p) 


3D(xq,p) 


(w,T^U^  »n)ds 


(7.4) 


This  equation  can  be  rewritten 


/  [(u,T^w^»n)  -  (w,T^u^  »n)  ]ds 

3D 


(7.5) 


=  /  [-(u,T^w^»n)  +  (w,T^-uJ  »n)  ]ds  . 

r 


The  functions  ♦  associated  to  the  solution  u  can  be 

written  in  the  neighborhood  of  z0: 


♦  (z)  =  aQ  +  a^ ( z-z0)  +  0((z-z0)2) 


(7.6) 
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we  get 


o 

& 

II 

+  b^ (z-z 0)  +  0(( z-z0)2) 

(7.7) 

<p(z)  = 

C(z)  -  z04' (z) . 

(7.8) 

.  3 )  and 

(7.6)-(7.8)  in  (7.5)  and 

letting  p  ♦  0 

•n)ds  - 

/  ( w,T^u^  »n)ds 

3D 

+  A)  + 

b^(<A  +  A)  +  a^(1  +  ie)B  +  a^ 

(1+k)BJ.  (7.9) 

of  (3. 

4)— (3-6)  we  get 

T11(£o} 

=  Re(a,|  +  a^  -  b^  ) 

(7.10) 

t22(x0) 

=  Re(a^  +  a^  +  b^ ) 

(7.11  ) 

t1  2  (~q) 

=  In  b|  . 

(7.12) 

By  proper  selection  of  A,  B  we  can  obtain  that  the  right 
hand  aide  of  (7-9)  be  T1>;j*  Note  thab  any  choice  of 

and  £*  in  (7.1)  and  (7.2)  does  not  change  the  right  hand 
side  of  (7*9). 

In  our  problem  when  the  tractions  are  prescribed  at  8T),  the 
function  g(x)  =  T^*n  is  given.  (7.9)  can  therefore  be 
written  in  the  form 


P  = 


/  ( u« ,  T  ^  •  n )  d  s  -  /  (w,g)ds 

3D  3D 


(7.13) 


where  P  is  (for  proper  choice  of  A,  B)  the  exact  value  of  the 
stress  component  at  x  =  Xg.  Of  course  Uq  is  not  known  but 
u-cm  is.  Therefore  we  define 


/  (u„_,fT^w^  *n)ds  -  /  (w,g)ds 


(7.14) 


By  subtracting  (7.13)  (7.14),  the  error  in  the  extracted 
functional  Ppp  (provided  that  integrals  are  evaluated 
exactly)  is 

P  -  PPE  3  /  (*0  “  HPE’T[wj-n)ds* 


(7.15) 


Let  us  analyze  now  (7.15).  To  this  end,  let  v  =  (v1 ,v2)  be 

the  (exact)  solution  of  the  problem  when  tractions  T^w^»n  are 
prescribed  at  3D.  v  4%  because  y  is  singular  at  x  =  x^, 

but  y  is  not.  Existence  of  v  is  guaranteed  because 
satisfy  the  equilibrium  condition.  We  can  write 

J  {x±0  *  HpE)'TW,n)ds  =  2W(uq  -  Upjj.v)  (7.16) 

O  U 

where  W(u,v)  is  the  usual  energy  scalar  product  associated 
with  W(u)  defined  in  (3.1).  Using  one  of  the  basic  property 
of  the  finite  element  method,  namely 


V(U  “  UpE.ipE)  =  °* 

(7.17) 

we  obtain  from 

(7.15)  (7.16) 

and  hence 

F  -  FpE  *  2W(u0  -  Upg ,  y  -  Vpg) 

-  fpe' 

1  <  2iuQ  -  Upgigiv  -  ypEig  . 

(7.18) 

So  far  we  did  not  discuss  the  choice  of  $#(z)  and  £*(z). 
(7.18)  shows  that  and  should  be  selected  so  that 

ly  -  vpE»  is  at  least  of  the  order  of  iuQ  -  upEi . 

If  “  Ype’e  S  Cl-  "  Hpe^  we  get  *P  “  PPE*  * 

CIUq  -  Upglg  <  C(W(uQ)  -  W(upg))  and  the  rate  of  convergence 

is  twice  that  of  the  rate  of  the  error  measured  in  the  energy 
norm.  Note  that  inequality  (7.18)  is  upper  bound  which 
neglects  possible  cancellation  in  the  energy  integral. 

8.  SELECTION  OP  THE  EXTRACTION  FUNCTION 

When  Xq  is  not  close  to  the  boundary  of  D,  then  we  can 

select  =  0.  When  Xq  is  close  to  3D,  then 

and  should  be  selected  so  that  T^w^»n  =0  on  that 

part  of  the  boundary  which  is  close  to  Xq.  Otherwise,  we 
would  not  achieve  that  iv  -  Vp^ig  will  be  small. 

In  the  following  we  outline  briefly  the  procedure  for 

constructing  and  so  that  T^w^»n  *  0  on  the  crack 

surfaces.  To  simplify  the  notation  we  will  write  f  instead 
of  +  ,  etc. 

Define  an  auxiliary  function  a(z)  on  D 
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q(z)  =  *(z)  +  z*’(z)  +  ?(z)  (8.1) 

Using  (3.4)  and  (3*5)  the  tractions  on  the  crack  surface  can  be 

written  as  follows 

t22(z+)  -  ir12(z+)  =  #(z+)  +  a( z_)  (8.2a) 

x22(z_)  -  iT1 2(z_)  =  ♦ ( z_ )  +  n(z+)  (8.2b) 

where  z+  and  z_  respectively  denote  the  upper  and  lower 
surface  of  the  crack.  Using  (7.1  )  —  (7 . 3 )  we  get 

n(z)  =  -X(z-z0)"2+2X(z-z0)(z-z0)“^-B(z-z0)~2+n#(z)  (8.3) 

where 

Q#(z)  a  ¥#(z)  +  z#;(z)  +  <»#(z).  (8.4) 

Setting 

fl#(z)  =  **(z)  (8.5) 

we  obtain 

t22(z+)  -  iT-j2^z+^  =  ♦  **(z+)  +  (9.6a) 


where 

Q ( z)  =  -A(z-zq)_2  +  X(z-z0)"2 

"  (8.7) 

+  2X(z  -  zQ ) ( z  -  z0)"3  -  B(z  -  z0)-2  . 

Note  that  Q(z+)  =  Q(z_).  Similarly 

< 

t22(z_)  -  ix12(z_)  =  Q(z+)  +  **(z_)  +  *#(z+)*  (8.6b) 

Now  we  select  •»  so  that 

♦#(z+)  +  e#(z_)  =  -Q(z+).  (8.8) 

(8.5)  and  (8.8)  define  now  4*  and  **.  By  this  selection  we 
achieve  that  t22(z+)  =  t12(z+)  =  T22^z-^  *  T12^Z+^  *  °*  Th® 

relation  (8.8)  can  be  easily  achieved.  For  example,  for 
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we  get 


4  ( z^2  +  z^2  ) 2 


*(z+)  +  *(z_) 


1 


which  is  one  terra  in  (8.7).  Consequently  we  get  the  other 
terras  and  combining  them  (8.8)  is  achieved. 


9.  NUMERICAL  PERFORMANCE  OF  THE  EXTRACTION  TECHNIQUE 

We  now  present  the  results  of  computational  experiments  based 
on  our  model  problem  and  the  extraction  function  described  in 
Section  8  (using  $*,  £#). 

Fig.  9.  shows  the  results  analogous  to  those  shown  in  Fig.  5 
but  stress  components  was  computed  by  the  extraction 

technique.  The  slope  shown  in  the  figure  shows  the  rate  u  =  1 
(i.e.,  the  rate  of  the  convergence  of  the  energy  and  not  the 

energy  norm).  For  comparison  the  error  e^2  for  mesh  2 

computed  directly  (see  Fig.  5b)  is  shown  also  in  Fig.  9.  Fig. 
10  shows  the  isometric  drawings  (in  the  same  scale  as  in 
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Figure  9 

The  relative  error  of  computed  by  postprocessing  of  the 

p-version  for  Mesh  1  and  Mesh  2.  (1)  e^  ,  (2)  ej^t  (3)  e^2 


The  behavior  of  the  error  x 22  comPu-ted  by 

postprocessing  of  the  p-version. 


computed  by  the  postprocessing 


Pig.  6)  of  the  error  in  t2_ 
technique. 

p 

Table  7  shows  the  relative  error  e^  in  the  stresses  at 

the  point  (.2^,.25)  computed  by  theipostprocessing  technique 
taking  i#  =  e#  =  0  (because  the  point  is  not  close  to  the 
boundary;.  This  data  should  be  compared  with  the  results  of 
Table  5  and  6. 


TABLE  7 

The  relative  error  efj  in  the  stresses  Tji  at  the  point 
(.25,. 25)  computed  by  postprocessing.  J 


No.  of 

elements 

N 

eR 

e11 

eR 

e22 

eR 

e12 

16 

51 

20.91* 

22.51* 

12.70* 

0) 

> 

43 

143 

11.40* 

13.07* 

8.29* 

£■§ 

106 

221 

4.60* 

5.92* 

3.74* 

319 

617 

1.47* 

2.01* 

1.31* 

5 

16 

51 

20.91* 

22.51* 

12.70* 

e 

64 

167 

12.51* 

14.88* 

9-86* 

.o  ,c 

256 

591 

6.87* 

8.90* 

6.28* 

■p 

We  see  that  for  adaptive  meshes  the  error  is  of  order  W1  and 

for  uniform  meshes  of  order  NJ/2  .  Similarly,  as  in  the  case  of 
the  p-version  we  see  an  orderly  convergence  with  the  rate  as 
the  square  of  the  error  measured  in  the  energy  norm  (as 
theoretically  expected). 


10.  CONCLUSIONS 

The  shown  computations  are  characteristic  in  the  following 
way.  The  convergence  in  the  energy  norm  is  monotonic  and  very 
orderly.  For  the  smooth  solution  the  p-version  is  especially 
effective.  For  unsmooth  solutions  the  refinement  of  the  meshes 
in  the  h-version  is  very  essential. 

The  convergence  of  stresses  in  a  fixed  point  is  very  "chaotic," 
the  accuracy  in  various  components  can  be  very  different.  The 
rate  of  convergence  of  the  postprocessed  values  are  as  the 
square  of  the  error  measured  in  the  energy  norm.  In  the  case 
of  the  h-version,  uniform  (or  piecewise  uniform)  meshes  and 
smooth  solution  the  supsreonvergence  occurs  in  the  center  of 

the  elements.  The  rate  is  h^  log  h,  i.e.,  effectively  as  the 
square  of  the  error  in  energy  norm  (h).  Therefore,  the  gain 


for  the  elements  of  degree  1  is  not  in  the  rate  of 
convergence  of  the  postprocessed  value  but  is  in  the 
magnitude.  (For  p  ">  1  the  gain  of  the  postprocessing  appears 
also  in  the  rate.) 

The  postprocessing  is  especially  important  for  the  p-version, 
although  it  is  also  essential  for  the  h-version  especially  for 
unsmooth  solutions  and  for  general  meshes. 

11.  EFFECTIVITY  OF  THE  POSTPROCESSING  TECHNIQUE 

In  the  introduction  we  raised  a  number  of  questions  concerning 
the  postprocessing.  We  now  briefly  address  these  question  in 
the  light  of  our  results.  Detailed  analysis  will  be  made  in  a 
forthcoming  paper. 

1 )  It  is  cost  effective  not  to  save  computational  effort 
on  a  postprocessing  procedure  especially  when  not  an  excessive 
number  of  extractions  is  made.  The  cost  of  obtaining  reliable 
and  accurate  values  by  postprocessing  is  much  smaller  than  to 
obtain  comparable  accuracy  by  increasing  p  in  the  p-version 
or  refine  the  meshes  in  the  h-version.  The  postprocessing 
usually  removes  very  reliably  the  "chaotic"  behaviour  of  the 
errors  in  stresses.  The  effectivity  of  the  postprocessing  is 
characterized  by  higher  rate  of  convergence  than  in  the  energy 
norm. 

2)  The  rate  of  convergence  as  the  square  of  the  rate  of 
the  error  in  the  energy  norm  is  theoretically  the  maximal  one 
which  can  be  directly  extracted.  The  postprocessing  technique 
we  outlined  leads  to  this  rate. 

3)  Developoment  and  implementation  of  the  postprocessing 
techniques  in  finite  element  programs  is  practically  not  a  very 
simple  task.  We  mention  some  aspects: 

a)  A  number  of  extraction  functions  must  be 
developed.  Although  many  analytical  solutions  of  special 
problems  are  very  helpful  for  such  development,  the 
general  approach  especially  for  nonhoraogeneous  material 
still  needs  further  research. 

b)  Special  care  must  be  excercised  in  the  numerical 
evaluation  of  integrals  because  the  extraction  function 
can  have  singular  character. 

c)  The  postprocessing  technique  for  nonlinear 
problems  could  be  especially  important  but  additional 
research  is  necessary. 
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